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Abstract 

Hysteresis  exhibited  by  smart  materials  hinders  their  wider  applicability  in  actuators  and  sensors.  In  this 
paper  methods  are  studied  for  recursive  identification  and  adaptive  inverse  control  of  smart  material  actuators, 
where  a  Preisach  operator  with  a  piecewise  uniform  density  function  is  used  to  model  the  hysteresis.  Persistent 
excitation  conditions  for  parameter  convergence  are  discussed  in  terms  of  the  input  to  the  Preisach  operator. 

Two  classes  of  recursive  identification  schemes  are  explored,  one  based  on  the  hysteresis  output,  the  other 
based  on  the  time  difference  of  the  output.  Asymptotic  tracking  for  the  adaptive  inverse  control  method  is 
proved,  and  the  condition  for  parameter  convergence  is  given  in  terms  of  the  reference  trajectory.  Practical 
implementation  issues  are  also  investigated.  Simulation  and  experimental  results  based  on  a  magnetostrictive 
actuator  are  used  to  illustrate  the  approach. 

1  Introduction 

Smart  materials,  e.g.,  magnet ostrictives,  piezoelectrics,  and  shape  memory  alloys  (SMA),  display  coupling  of 
elastomechanics  with  electromagnetic/thermal  influences.  Hence  these  materials  have  built-in  sensing/actuation 
mechanisms.  Hysteresis  widely  existing  in  smart  materials,  however,  makes  the  effective  use  of  smart  material 
actuators  and  sensors  quite  challenging.  Control  of  hysteresis  in  smart  materials  has  attracted  attention  in  recent 
years  [1].  Hysteresis  models  can  be  roughly  classified  into  physics-based  models  and  phenomenological  models. 
The  most  popular  phenomenological  hysteresis  model  used  for  smart  materials  has  been  the  Preisach  model 
[2,  3,  4,  5,  6,  7,  8].  A  similar  type  of  operator  called  Krasnoserskii-Pokrovskii  (KP)  operator  has  also  been  used 
[9,  10]. 

Inverse  compensation  is  a  fundamental  approach  in  coping  with  hysteresis,  where  one  aims  to  cancel  out  the 
hysteresis  effect  by  constructing  a  right  inverse  of  the  hysteresis  operator  [3,  11,  12,  10,  7,  8].  The  performance  of 
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inverse  compensation  is  susceptible  to  model  uncertainties  and  to  errors  introduced  by  (inexact)  inverse  algorithms. 
Hysteretic  behaviors  of  smart  materials  often  vary  with  time,  temperature  and  some  other  ambient  conditions.  To 
combat  this  problem,  a  robust  control  framework  was  proposed  by  combining  inverse  compensation  with  l\  control 
theory  in  [13].  An  alternative  approach  is  adaptive  inverse  control  [11,  14,  15].  Tao  and  Kokotovic  developed 
an  adaptive  inverse  control  scheme  for  a  class  of  hysteresis  models  with  piecewise  linear  characteristics  [11].  An 
adaptive  identification  and  inverse  compensation  method  was  studied  for  the  KP  operator  and  applied  to  a  SMA 
actuator  by  Webb  et  al.  [14].  In  [15]  an  adaptive  inverse  scheme  was  presented  for  piezoelectric  actuators,  where 
the  Prandtl-Ishlinskii  hysteresis  operator  was  used  for  hysteresis  modeling. 

This  paper  is  devoted  to  the  study  of  recursive  identification  and  adaptive  inverse  of  hysteresis  in  smart 
materials.  The  following  highlights  the  major  differences  of  this  work  from  the  previous  works  reported  in  [11,  14, 
15]:  (1)  the  Preisach  operator  is  used  as  the  hysteresis  model;  (2)  the  persistent  excitation  (P.E.)  conditions  for 
the  convergence  of  parameter  identification  are  studied  in  terms  of  the  input  to  the  hysteresis  operator;  and  (3) 
the  asymptotic  tracking  property  of  the  adaptive  inverse  control  algorithm  is  proved,  and  the  issue  of  parameter 
convergence  is  discussed  in  terms  of  the  reference  trajectory.  Experimental  results  based  on  a  magnetostrictive 
actuator,  together  with  extensive  simulation  results,  are  used  to  examine  the  effectiveness  of  the  proposed  approach. 

The  remainder  of  the  paper  is  organized  as  follows.  The  Preisach  operator  is  briefly  reviewed  in  Section  2. 
Recursive  identification  algorithms  are  studied  in  Section  3.  Two  classes  of  schemes  are  compared,  one  based  on 
the  output  of  the  hysteresis  model  and  the  other  based  on  the  time  difference  of  the  output.  The  adaptive  inverse 
control  scheme  is  discussed  in  Section  4.  Finally  some  conclusions  are  provided  in  Section  5. 


2  The  Preisach  Operator 


The  Preisach  operator  is  briefly  reviewed  in  this  section.  A  more  detailed  treatment  can  be  found  in  [16,  17,  18]. 
For  a  pair  of  thresholds  (/?,  a)  with  (3  <  a,  consider  a  simple  hysteretic  element  7/3, ab  *]>  as  illustrated  in  Fig.  1. 
Let  C([0,T])  denote  the  space  of  continuous  functions  on  [0,T].  For  u  G  C([0,T])  and  an  initial  configuration 
C  G  {  —  1,1},  l 0  =  7 /3,a[rt,C]  is  a  function  from  [0,  T]  to  {  —  1, 1}  defined  as  follows  [17]: 


"(0)  =  < 


if  u( 0)  <  (3 
if  (3  <  u{ 0)  <  a 
if  i/(0)  >  a 


and  for  t  G  (0,T],  setting  Xt  =  {r  G  (0,  t]  :  u(r)  =  (3  or  ct}, 


"(f)  =  < 


"(0) 

-1 

1 


if  Xt  =  0 

if  Xt  0  and  u( maxlt)  =  (3 
if  Xt  7^  0  and  u{ max  At)  =  a 


This  operator  is  sometimes  referred  to  as  an  elementary  Preisach  hysteron  (it  is  called  hysteron  hereafter  in 
this  paper),  since  it  is  a  building  block  of  the  Preisach  operator.  Define  Vo  =  {(/?,  a)  G  M2  :  f3  <  a}.  Vo  is 
called  the  Preisach  plane ,  and  each  (/?,  a)  G  Vo  is  identified  with  the  hysteron  7/3, a-  For  u  G  C([0,  T])  and  a  Borel 
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Fig.  1:  An  elementary  Preisach  hysteron  7/3, ab  *]• 

measurable  initial  configuration  Co  of  all  hysterons,  Co  :  Vo  — >  {  —  1,1},  the  output  of  the  Preisach  operator  T  is 
defined  as  [17]: 

y(t)  =  T[u,(o](t)  =  [  Ii((3,a)7i3,a[u,(o03,a)\(t)d(3da,  (1) 

Jv o 

where  the  weighting  function  /a  is  often  referred  to  as  the  Preisach  function  [16]  or  the  density  function  [18]. 
Throughout  the  paper  it  is  assumed  that  fi  >  0.  Furthermore,  to  simplify  the  discussion,  assume  that  /a  has  a 
compact  support,  i.e. ,  fi(/3,a)  =  0  if  f3  <  /3o  or  a  >  ao  for  some  {3o,olo-  In  this  case  it  suffices  to  consider  a  finite 

triangular  area  in  the  Preisach  plane  V  =  {(/?,  a)  G  Vo\(3  >  An  <  <^o}  (see  Fig.  2(a)). 

The  memory  effect  of  the  Preisach  operator  can  be  captured  by  the  memory  curves  in  V.  At  time  £,  V  can  be 
divided  into  two  regions: 

V+(t)  =  {(/?,  a)  G  V\  output  of  73, a  at  t  is  +  1}, 

V-(t)  =  { (/?,£*)  G  V\  output  of  jpjOL  at  t  is  —  1}. 

Now  assume  that  at  some  initial  time  to,  the  input  u(to)  =  uo  <  An  Then  the  output  of  every  hysteron  is  —  1. 
Therefore  V-(to)  =  V,  V+(to)  =  0  and  this  corresponds  to  the  “negative  saturation”  (Fig.  2(b)).  Next  assume 
that  the  input  is  monotonically  increased  to  some  maximum  value  at  t\  with  u(t\)  =  u\.  The  output  of  73, a  is 
switched  to  +1  as  the  input  u(t)  increases  past  a.  Thus  at  time  £i,  the  boundary  between  V~(t\)  and  V+(t\)  is  the 
horizontal  line  a  =  u\  (Fig.  2(c)).  Next  assume  that  the  input  starts  to  decrease  monotonically  until  it  stops  at  1 2 
with  ufo)  =  U2-  It’s  easy  to  see  that  the  output  of  7^  becomes  —1  as  u(t)  sweeps  past  /?,  and  correspondingly, 
a  vertical  line  segment  (3  =  U2  is  generated  as  part  of  the  boundary  (Fig.  2(d)).  Further  input  reversals  generate 
additional  horizontal  or  vertical  boundary  segments. 

From  the  above  illustration,  each  of  V-  and  T+  is  a  connected  set  [16],  and  the  output  of  T  is  determined  by 
the  boundary  between  V-  and  T+  if  there  is  no  Preisach  weighting  mass  (or  impulse  in  fi)  on  the  curve.  The 
boundary  is  called  the  memory  curve.  The  memory  curve  has  a  staircase  structure  and  its  intersection  with  the 
line  a  =  /3  gives  the  current  input  value.  The  set  of  all  memory  curves  is  denoted  T.  The  memory  curve  ^ 0  at 
t  =  0  is  called  the  initial  memory  curve  and  it  represents  the  initial  condition  of  the  Preisach  operator.  Hereafter 
the  initial  memory  curve  will  be  put  as  the  second  argument  of  the  Preisach  operator. 

Rate-independence  is  one  of  the  fundamental  properties  of  the  Preisach  operator: 

Theorem  2.1  (Rate-independence)  [17]  If  <f  :  [0,T]  — >  [0,T]  is  an  increasing  continuous  function  satisfying 
(f)(0)  =  0  and  (f)(T)  =  T,  then  for  u  G  C([0,T]);  r[wo^^0](t)  =  T[u,  x[>o\(4>(t)) ,  Vt  G  [0,T],  where  “o”  denotes 
composition  of  functions. 
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Fig.  2:  Memory  curves  in  the  Preisach  plane. 

3  Recursive  Identification  of  Hysteresis 

For  the  Preisach  operator,  the  “parameter”  is  the  Preisach  density  function.  A  classical  method  to  identify  the 
density  function  is  using  the  so  called  first  order  reversal  curves ,  detailed  in  Mayergoyz  [16].  A  first  order  reversal 
curve  can  be  generated  by  first  bringing  the  input  to  /?o,  followed  by  a  monotonic  increase  to  <a,  then  a  monotonic 
decrease  to  (3.  The  term  “first  order  reversal”  comes  from  that  each  of  these  curves  is  formed  after  the  first 
reversal  of  the  input.  Denote  the  output  value  as  /(/3,  a)  when  the  input  reaches  ft.  Then  the  density  fi(/3,a) 
can  be  obtained  as  /. i(/3,a )  =  \  •  Since  it  involves  twice  differentiation,  a  smooth  approximating  surface 

is  fit  to  the  data  points  in  practice  [3,  6].  As  pointed  out  in  [6],  deriving  the  density  by  differentiating  a  fitted 
surface  is  inherently  imprecise,  since  different  types  of  approximating  functions  lead  to  quite  different  density 
distributions.  A  second  approach  is  to  devise  the  input  sequence  carefully  and  derive  the  Preisach  weighting 
masses  (on  a  discretization  grid)  directly  from  the  output  measurements  [19].  This  scheme  is  very  sensitive  to 
measurement  noises  as  one  can  easily  see.  Hence  a  third  approach  is  to  identify  the  Preisach  weighting  masses  or 
the  weights  for  basis  functions  based  on  the  least  squares  method  [9,  10,  7].  All  of  the  aforementioned  methods 
are  used  off-line.  However,  it  has  been  observed  that  the  hysteretic  behaviors  of  smart  materials  are  dependent 
on  temperature  and  some  other  ambient  conditions,  and  often  vary  with  time  slowly.  Hence  it  is  desirable  to  have 
recursive,  on-line  schemes  for  the  identification  of  hysteresis  parameters. 


3.1  Discretization  of  the  Preisach  operator 

In  practice  the  Preisach  operator  needs  to  be  discretized  in  one  way  or  another  during  the  identification  process. 
A  natural  way  to  approximate  a  Preisach  operator  is  to  assume  that  inside  each  cell  of  the  discretized  Preisach 
plane,  the  Preisach  density  function  is  constant.  This  approximation  has  nice  convergence  (to  the  true  Preisach 
operator)  properties  under  mild  assumptions  [19]. 

Let  [umin,  Umax]  be  the  practical  input  range  to  the  hysteresis  operator,  which  is  often  a  strict  subset  of  [ao,(3o\. 
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For  the  hysteretic  behavior  one  can  focus  on  the  triangle  bounded  by  umin  and  umax  in  the  Preisach  plane,  since  the 
contribution  to  the  output  from  hysterons  outside  this  triangle  is  constant  [7].  Discretize  uniformly 

into  L  +  1  levels  (  called  discretization  of  level  L  in  this  paper),  where  the  discrete  input  levels  Ui,  l<i<L  +  l, 
are  defined  as 

U{  —  Umin  H-  if  l)An, 

where  An  =  .  The  cells  in  the  discretization  grid  are  labeled,  as  illustrated  in  Fig  3(a)  for  the  case  of 

L  =  4. 


a 


a 


(a) 


W5 


(b) 


Fig.  3:  Illustration  of  the  discretization  scheme  (L  =  4):  (a)  Labeling  of  the  disretization  cells;  (b)  Weighting 
masses  sitting  at  the  centers  of  cells. 


We  note  that  the  Preisach  operator  with  piecewise  uniform  density  is  still  an  infinite-dimensional  operator.  If 
one  assumes  that  the  Preisach  weighting  function  inside  each  cell  is  concentrated  at  the  center  as  a  weighting  mass 
(Fig.  3(b)),  the  corresponding  Preisach  operator  becomes  a  weighted  combination  of  a  finite  number  of  hysterons. 


3.2  Recursive  identification  schemes 

In  the  interest  of  digital  control,  the  discrete-time  setting  is  considered  in  this  paper.  One  way  to  obtain  a 
piecewise  uniform  density  is  to  first  identify  the  discrete  weighting  masses  (Fig.  3(b)),  and  then  distribute  each 
mass  uniformly  over  the  corresponding  cell  [13].  A  Preisach  operator  with  discrete  weighting  masses  is  easier  to 
analyze  than  a  Preisach  operator  with  a  piecewise  uniform  weighting  density;  however,  these  two  types  of  operators 
bear  much  similarity  and  essential  results  for  one  can  be  easily  translated  into  those  for  the  other.  Hence  recursive 
identification  of  Preisach  weighting  masses  is  first  studied,  and  then  the  extension  needed  for  identifying  the  density 
directly  is  briefly  discussed. 

In  this  paper  two  classes  of  identification  algorithms  are  examined,  one  based  on  the  hysteresis  output,  and 
the  other  based  on  the  time  difference  of  the  output  (called  difference-based  hereafter). 
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Output-based  identification 


The  output  y[n\  of  the  discretized  Preisach  operator  (corresponding  to  the  case  illustrated  in  Fig.  3(b))  at  time 
instant  n  can  be  expressed  as 

L  i 


y[n\  =  EE^'MTp 

i= 1  3= 1 


(2) 


where  Wij  [n]  denotes  the  state  (1  or  —1)  of  the  hysteron  in  cell  (i,j)  at  time  n,  and  u*-  denotes  the  hys- 
teron’s  Preisach  weighting  mass.  Stacking  Wij  [n]  and  v*-  into  two  vectors,  W[n\  =  [Wi[n]  •  •  •  Wx[n]]T  and 
is*  =  Wi  ’  ’  ’  vk\Ti  where  K  =  L^+1^  is  the  number  of  cells,  one  rewrites  (2)  as 


K 


y[n\  =  Y,Wk[nVk  =  W[n\Tv* 


k= 1 


(3) 


Let  ts[n\  =  [z>i  [n]  •  •  •  z>^[n]]T  be  the  estimate  of  z/*  at  time  n,  and  let 


K 

y\n\  =  y>»[n]fr[n]  =  W[n}Tv[n ] 

fe= l 


be  the  predicted  output  based  on  the  parameter  estimate  at  time  n.  The  gradient  algorithm  [20]  to  update  the 
estimate  is 


z>[n  +  1]  =  0[n\  —  7 


(j/M  -  y[n])W[n] 
W[n]TW[n] 


(4) 


where  0  <  7  <  2  is  the  adaptation  constant.  To  ensure  that  the  weighting  masses  are  nonnegative,  we  let 
[n  +  1]  =0  if  the  k- th  component  of  the  right  hand  side  of  (4)  is  negative.  Since  this  parameter  projection  step 
brings  the  parameter  estimate  closer  to  the  true  values,  it  does  not  invalidate  the  convergence  result  should  this 
step  were  absent  [20].  Hence  this  step  will  not  be  explicitly  considered  for  convergence  analysis  in  this  paper  (in 
particular,  in  the  proof  of  Theorem  4.1).  Define  the  parameter  error  D[n \  =  0[n\  —  is* .  Then 


z>[n  +  1]  =  F[n]z>[n], 


where  F[n\  =  Ik 


W[n]W[n]T 
W[n]TW[n]  ’ 


and  I k  represents  the  identity  matrix  of  dimension  K. 


(5) 


It  is  well-known  [20]  that  the  convergence  of  the  algorithm  (4)  depends  on  the  persistent  excitation  (P.E.) 
condition  of  the  sequence  W[n\.  The  sequence  W[n\  is  persistently  exciting  if,  there  exist  an  integer  N  >  0  and 
c[  >  0,  cf2  >  0,  such  that  for  any  no, 


n0+N— 1 

c[Ik  <  ^2 

n=rio 


W[n]W[n]T 

W[n]TW[n\ 


<  c'2Ik- 


(6) 


Due  to  the  equivalence  of  uniform  complete  observability  under  feedback  [20,  21],  from  (6),  there  exist  c\  >  0,  C2  >  0 
such  that  for  any  no, 

ci Ik  <  Gat(^o)  <  c2Ik:  (7) 


where  Gat(^o)  is  the  observability  grammian  of  the  system  (5)  defined  as 


Gat(^o) 


n0+W— 1 

E 

n=no 


4>[n,  no]TFF[n]lF[n]T4>[n,  no] 
W[n]TW[n] 
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and  $[71,  no]  is  the  state  transition  matrix  <L[n,  no]  =  Ylk=n0  F[k\-  It  can  be  shown  [20]  that  when  (7)  is  satisfied, 


\\u[n  +  N}\\  <  v'l  -  -  n  | /'/[;;] |  • 


from  which  exponential  convergence  to  is*  can  be  concluded. 


(8) 


Difference-based  identification 


An  alternate  way  to  identify  is*  is  using  the  time  difference  z[n\  of  the  output  y[n\,  where 

z[n\  =  y[n }  —  y[n  —  1]  =  (W[n\  —  W[n  —  1])T  v* . 


Let  y[n  }  and  y[n  —  1]  be  the  output  predictions  at  time  n  and  n  —  1  based  on  0[n  —  1],  respectively,  i.e. , 


y[n  ]  =  W[n]T0[n  —  1],  y[n  —  1]  =  W[n  —  1]T i >[n  —  1]. 

Define  z[n]  =  y[n~ ]  —  y[n  —  1]  =  (W[n]  —  W[n  —  l])Tz>[n  —  1].  Let  V[n\  be  the  time  difference  of  hysteron  states, 
V[n\  =  W[n]  —  W[n  —  1].  Then  we  can  obtain  the  following  identification  scheme  based  on  z\n\: 


v[n  +  1] 


ifF[n]^0 

v[n\  if  V[n\  =  0 


(9) 


As  in  the  output-based  scheme,  an  additional  parameter  projection  step  will  be  applied  if  any  component  of 
z>[n  +  1]  is  negative.  Similarly  one  can  write  down  the  error  dynamics  equation,  the  P.E.  condition  on  V[n\,  and 
the  convergence  rate  estimate  corresponding  to  the  difference-based  scheme  (9).  The  sequences  V[n\  and  W[n\ 
are  almost  equivalent  in  the  sense  that,  for  any  N  >  0,  {V[n]}^=1  can  be  constructed  from  {W[n\}^=  0,  and 
conversely,  {W\n]}^=1  can  be  constructed  from  W[ 0]  and  {V[n]}^=1.  However,  there  are  motivations  to  introduce 
the  difference-based  scheme  (9).  While  W[n\  has  components  ±1,  the  components  of  V[n\  are  ±2  or  0.  Often 
times  most  components  of  V[n\  are  0  since  V/Jn]  7^  0  only  if  the  k- th  hysteron  changed  its  state  at  time  n.  This 
has  two  consequences:  (1)  The  P.E.  condition  of  V[n\  is  easier  to  analyze  than  that  of  W[n\;  (2)  The  convergence 
of  the  difference-based  scheme  (assuming  that  P.E.  is  satisfied)  is  faster  than  that  of  the  output-based  scheme 
since  z[n\  carries  more  specific  information  about  v*. 

The  P.E.  condition  for  the  difference-based  algorithm  is  equivalent  to  that  spans  RK  since  V[n\ 

can  take  only  a  finite  number  of  possible  values.  It  is  of  practical  interest  to  express  the  P.E.  conditions  in  terms 
of  the  input  u[n\  to  the  hysteresis  operator  (in  which  case  u[n\  is  said  to  be  P.E.).  Recall  that  u[n\  takes  values  in 
a  finite  set  { U{ ,  1  <  i  <  L  +  l}.  To  avoid  ambiguity  one  should  understand  that  the  input  to  the  Preisach  operator 
is  monotonically  changed  from  u[n—  1]  to  u[n\.  In  the  analysis  below  it  is  assumed  that  the  input  does  not  change 
more  than  one  level  during  one  sampling  time.  The  assumption  is  not  restrictive  considering  the  rate-independence 
of  the  Preisach  operator,  but  it  helps  to  ease  the  presentation. 


Theorem  3.1  (Necessary  condition  for  P.E.)  If  {V[n]}  is  P.E.,  then  there  exists  N  >  0,  such  that  for  any 
no,  for  any  i  G  {1,  2,  *  •  *  ,L};  u[n\  achieves  a  local  maximum  at  or  a  local  minimum  at  U{  during  the  time 
period  [no,  no  +  N  —  1] . 
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Proof.  Let’s  call  a  hysteron  active  at  n  if  it  changes  state  at  time  n.  Since  the  input  changes  at  most  one  level 
each  time,  if  u[n]  >  u[n  —  1],  the  set  of  active  hysterons  must  have  the  form  =  {(i,  j),  (i,  j  +  1),  •  •  •  ,  (i,  i)}  for 
some  i,  j  with  1  <  i  <  L  and  1  <  j  <  i  (refer  to  the  labeling  scheme  in  Fig.  3(a)),  and  the  components  of  V[n\ 
corresponding  to  elements  of  S^-  are  2  and  other  components  equal  0.  Similarly,  if  u[n\  <  u[n  —  1],  the  set  of  active 
hysterons  has  the  form  Sf-  =  {(j,  j),  (j  +  1,  j),  •  •  •  ,  (i,  j)}  for  some  i,j,  and  the  components  of  V[n\  corresponding 
to  elements  of  Sf-  are  —2  and  other  components  equal  0. 

If,  for  some  i' ,  Ui>+ \  is  not  a  local  maximum  and  u\  is  not  a  local  minimum,  or  Sff  it  will  not  become  the 
set  of  active  hysterons  during  [no,  no  # N  —  1].  Analysis  on  the  Preisach  plane  reveals  that  the  contribution  to  the 
output  from  the  hysteron  (i',if)  cannot  be  isolated,  and  hence  {VnYn=n^+ l  does  not  span  □ 


Remark  3.1  From  Theorem  3.1 ,  for  a  Preisach  operator  with  discretization  level  L,  it  is  necessary  that  the  input 
u[n\  has  L  reversals  at  different  input  levels  for  parameter  convergence.  This  is  in  analogy  to  ( but  remarkably 
different  from)  the  result  for  linear  systems,  where  the  input  is  required  to  have  at  least  n  frequency  components 
for  identification  of  n  parameters  [20,  21  ]. 


Theorem  3.1  implies  that  the  input  levels  u\  and  ul+ i  must  be  visited  for  P.E.  to  hold.  When  the  input  hits 

ui,  all  hysterons  have  output  —1  and  the  Preisach  operator  is  in  negative  saturation;  similarly,  when  the  input 

hits  ul+ i,  the  Preisach  operator  is  in  positive  saturation.  For  either  case  all  the  previous  memory  is  “erased” 

and  the  operator  is  “reset”.  Starting  from  these  reset  points,  one  can  keep  track  of  the  memory  curve  f>[n\  (the 

state  of  the  Preisach  operator)  according  to  the  input  u[-].  Consider  an  input  sequence  {u[n\}rfb=ria,  na  <  n^.  If 

there  exist  n\,  712,^3  and  714  with  na  <  n\  <  712  <  ns  <  714  <  n 5  such  that  the  memory  curve  V^i]  = 

and  f>[n2]  =  ^[714],  we  can  obtain  another  input  sequence  {u'[n]}7)f=na  by  swapping  the  section  {77[77]}^Lni  with 

P  E 

the  section  { u[n ]}™ln3.  We  write  {u[n]}rfb=na  =  {u'[n]\nb=na  (called  equivalent  in  terms  of  P.E.)  since  the 
two  sequences  carry  same  excitation  information  for  the  purpose  of  parameter  identification.  The  set  of  all  input 
sequences  obtained  from  }  as  explained  above  (with  possibly  zero  or  more  than  one  swappings)  form  the 

P.E.  equivalent  class  of  {u[n\ }™Lna}>  denoted  as  {^M}n=na}-  Note  that  in  particular,  {^M}nLna}  G  Mn]}n=na}- 
We  are  now  ready  to  present  a  sufficient  condition  for  P.E.  in  terms  of  the  input  u[n\. 


Theorem  3.2  (Sufficient  condition  for  P.E.)  If  there  exists  N  >  0,  such  that  for  any  no,  one  can  find 
{U[n]}n^n^  1  satisfying  the  following:  there  exist  time  indices  no  <  na  <  nx  <  mf  <  n2  < 
n2  <  •  •  •  <  nf  <  7i+  <  •  •  •  <  nb  <  n0  +  N  —  1  or  no  <  na  <  mf  <  nf  <  n2  <  nf  <  •  •  •  <  77+  <  rif  <  •  •  •  <  n^  < 
no  +  N  —  1,  such  that  uf[nf]  is  a  local  maximum  and  u'[nf]  is  a  local  minimum  of  {u'[n]}7ff=ria  for  each  i,  these 
local  maxima  and  minima  include  all  input  levels  { }^=i,  and  either 

(a)  {u'[n)~]}  is  non-increasing,  u'[n)~]  >  u'[n\  for  nf  <  n  <  n\>,  u'[n[~]  differs  from  id  [71^]  by  no  more  than  Au, 
and  {u'[nf}}  is  non- decreasing,  u'[nf]  <  u'[n\  for  nf<n<  n^,  u'[nf]  differs  from  ur[nf+1\  by  no  more  than  Au; 
or 

(b)  {u'[n)~]}  is  non- decreasing,  u'[n)~]  <  u'[n\  for  n+  <  n  <  n b,  u'[n)~]  differs  from  by  no  more  than  Au, 

and  [u'[nf]}  is  non-increasing,  u'[nf]  >  u'[n\  for  nf<n<  n \>,  u'[nf]  differs  from  u'[nf+f\  by  no  more  than  Au, 

then  V[n\  corresponding  to  u[n]  is  P.E. 


8 


Proof.  Construct  a  new  input  sequence  u[n\  based  on  u'[n\  which  achieves  the  local  maxima  {u'[n+]}  and  the 
local  minima  {u'[nf]}  with  the  same  order  as  in  uf[n],  but  u[n\  varies  monotonically  from  a  maximum  to  the  next 
minimum  or  from  a  minimum  to  the  next  maximum.  For  such  an  input,  memory  curve  analysis  on  the  Preisach 
plane  reveals  that  the  corresponding  {H[n]}  spans  ~RK .  From  the  way  u[n\  is  constructed  and  the  conditions  given 
in  the  theorem,  any  vector  in  {H[n]}  must  also  be  present  in  {^'[n]}  corresponding  to  u'[n\.  Hence  {H'[n]}  is  P.E. 
Finally  P.E.  of  {H[n]}  follows  since  {u'[n\ belongs  to  the  P.E.  equivalent  class  of  □ 

Theorem  3.2  is  not  conservative,  and  it  covers  a  wide  class  of  P.E.  inputs.  For  example,  it  can  be  easily 
verified  that  a  (periodic)  first  order  reversal  input  (see  Fig.  4(a)  for  case  L  =  4),  which  has  been  widely  used 
for  identification  of  Preisach  density  function  [16],  and  a  (periodic)  oscillating  input  with  decreasing  amplitude 
(Fig.  4(b)  for  case  L  =  4)  both  satisfy  the  conditions  in  Theorem  3.2,  and  are  thus  P.E.  In  these  two  cases,  u[n\ 
itself  satisfies  the  conditions  imposed  for  u'[n\  in  the  theorem.  Fig.  5  shows  an  example  where  one  can  conclude 
the  P.E.  of  a  periodic  u[n\  by  inspecting  a  P.E.  equivalent  input  u'[n\.  Note  that  Theorem  3.2  does  not  require 
u[n\  to  be  periodic,  although  periodic  examples  are  chosen  here  for  easy  illustration. 


(a)  (b) 


Fig.  4:  Examples  of  P.E.  inputs  (L  =  4,  showing  one  period):  (a)  The  first  order  reversal  input;  (b)  An  oscillating 
input  with  decaying  amplitude. 


Fig.  5:  An  example  of  P.E.  input  (L  =  4,  showing  one  period).  The  input  rd[n],  P.E.  equivalent  to  u[n\,  is  obtained 
by  swapping  two  sections  A  —  B  and  A'  —  B'  of  u[n\. 

3.3  Comparison  of  recursive  identification  schemes 

In  this  subsection  the  output-based  scheme  is  compared  with  the  difference-based  one  through  simulation.  As 
shown  in  (8),  the  minimum  eigenvalue  of  the  observability  grammian  (i.e.,  c\  in  (7))  is  directly  related  to  the 
convergence  rate  of  the  output-based  scheme.  The  same  statement  holds  for  the  difference-based  scheme  provided 
that  W[n\  is  replaced  with  V[n\  in  the  related  equations.  In  Table  1  we  list  the  corresponding  y/1  —  ~c[  (the  bound 
on  parameter  error  drop  over  one  period)  under  the  two  gradient  schemes  (with  7  =  1)  for  different  discretization 
levels  L  with  the  (periodic)  first  order  reversal  input.  From  Table  1,  the  difference-based  scheme  converges  faster 
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as  expected.  Simulation  has  been  conducted  for  the  case  L  =  10.  Fig.  6(a)  compares  the  decrease  of  parameter 
error  over  periods  for  the  two  algorithms  when  there  is  no  measurement  noise,  and  the  conclusion  is  consistent 
with  Table  1. 

Table  1:  Comparison  of  convergence  rates  for  the  output-based  algorithm  and  the  difference-based  algorithm. 


L 

V 1  -  Ci 

(Output-based) 

VI -ci 

(Difference-based) 

5 

0.9631 

0.9399 

10 

0.9908 

0.9784 

15 

0.9958 

0.9874 

20 

0.9976 

0.9912 

25 

0.9985 

0.9933 

(a)  (b) 

Fig.  6:  Comparison  of  parameter  convergence  for  the  output  based  algorithm  and  the  difference-based  algorithm, 
(a)  Case  I:  noiseless  measurement;  (b)  Case  II:  noisy  measurement. 


Despite  the  apparent  advantage  of  faster  convergence,  the  difference-based  scheme  is  more  sensitive  to  the 
measurement  noise:  the  noise  gets  magnified  when  one  takes  the  output  difference  (analogous  to  taking  the 
derivative  of  a  noisy  continuous-time  signal),  and  the  disturbance  is  shared  only  among  the  active  hysterons. 
Simulation  in  Fig.  6(a)  is  re-conducted  where  a  noise  is  added  to  the  output,  the  noise  magnitude  being  4%  of  the 
saturation  output  of  the  Preisach  operator.  From  Fig.  6(b),  in  this  case,  the  parameter  error  will  not  converge  to 
zero  under  either  of  the  two  algorithms.  However,  the  ultimate  identification  error  of  the  output-based  algorithm 
is  much  lower  than  that  of  the  difference-based  scheme. 


3.4  Experimental  results  based  on  a  magnetostrictive  actuator 

Having  discussed  the  methods  for  recursive  identification  of  weighting  masses  for  a  Preisach  operator,  we  now 
point  out  how  to  change  the  previous  algorithms  to  identify  the  (piecewise  uniform)  Preisach  density  directly.  In 
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this  case  the  input  is  no  longer  limited  to  a  finite  set  of  values;  instead  it  can  take  any  value  in  [umin,  Umax].  Each 
component  Wij  [n]  of  W[n\  no  longer  represents  the  state  (1  or  —1)  of  the  hysteron  at  the  center  of  the  cell  (i,  j); 
instead  it  represents  the  signed  area  of  the  cell: 

Wij  [n]  =  area  of  C^-  [n\  —  area  of  CT-  [n\ , 

where  C±-  [n]  =  {(/ 6,  a )  G  cell  (i,j)  :  the  output  of  7£?Q,  at  time  n  is  db  1}.  The  definition  for  V[n\  remains 
the  same:  V[n\  =  W[n\  —  W[n  —  1].  The  P.E.  conditions  (Theorems  3.1  and  3.2)  presented  for  the  case  of 
weighting  masses  can  be  extended  to  the  case  of  piecewise  uniform  densities  in  a  straightforward  manner  with 
minor  modifications. 

Experiments  have  been  conducted  on  a  magnetostrictive  actuator  to  examine  the  identification  schemes.  Mag¬ 
netostriction  is  the  phenomenon  of  strong  coupling  between  magnetic  properties  and  mechanical  properties  of  some 
ferromagnetic  materials  (e.g.,  Terfenol-D):  strains  are  generated  in  response  to  an  applied  magnetic  field,  while 
conversely,  mechanical  stresses  in  the  materials  produce  measurable  changes  in  magnetization.  Magnetostrictive 
actuators  have  applications  to  micro-positioning,  robotics,  ultrasonics,  vibration  control,  etc.  Fig.  7  shows  a  sec¬ 
tional  view  of  a  Terfenol-D  actuator  manufactured  by  Etrema  Products,  Inc.  By  varying  the  current  in  the  coil, 
one  varies  the  magnetic  field  in  the  Terfenol-D  rod  and  thus  controls  the  motion  of  the  rod  head.  Like  other  smart 
material  actuators,  magnetostrictive  actuators  display  strong  input-output  hysteresis  [22]. 


Stainless  Steel  Push  Rod 


Fig.  7:  Sectional  view  of  a  Terfenol-D  actuator  [23] (Original  source:  Etrema  Products,  Inc.). 

When  operated  in  a  low  frequency  range  (typically  below  5  Hz),  the  magnetostrictive  hysteresis  is  rate- 
independent  and  can  be  modeled  by  (see  [7]) 

{H  =  Co  I  +  H\)ias 

M  =  T[H^  0]  •  (10) 

D  =  hod  As  M2 

Here  I  is  the  current  input  and  D  is  the  displacement  output.  The  magnetic  field  H  along  the  rod  direction  is 
related  to  /  linearly,  where  Co  >  0  is  the  coil  factor  and  Hbias  is  the  bias  field  necessary  for  generating  bi-directional 
strains.  The  bulk  magnetization  M  along  the  rod  direction  is  related  to  D  via  a  square  law,  and  the  constants  lrod, 
As,  Ms  denote  the  rod  length,  saturation  magnetostriction  and  saturation  magnetization,  respectively.  Hence  the 
magnetostrictive  hysteresis  is  essentially  captured  through  the  ferromagnetic  hysteresis  between  M  and  i7,  which 
is  modeled  by  the  Preisach  operator  T. 

In  our  experiments  the  actuator  displacement  is  measured  with  a  LVDT  sensor,  which  (after  low-pass  filtering) 
has  a  precision  of  ±0.5/im.  The  DSpace  ControlDesk  is  used  to  send  control  commands  and  collect  data.  The 
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following  parameters  are  available  from  the  manufacturer:  Co  =  1.54  x  104/m,  Ms  =  7.87  x  105A/m,  lrod  = 
5.13  x  10-2m,  and  the  following  parameters  can  be  identified  relatively  easily:  As  =  1.3  x  10-3,  =  1.12  x 

104A/m.  The  current  input  is  limited  to  the  range  [— 0.7A,  1.3A].  “Practical”  negative  saturation  Mmin  and 
positive  saturation  Mma;r  can  be  obtained  by  measuring  the  displacements  corresponding  to  I  =  —OTA  and  1.3A, 
respectively.  The  constant  contribution  to  the  Preisach  operator  (refer  to  the  discussion  in  Subsection  3.1  )  is 
evaluated  as  . 

A  periodic  first  order  reversal  current  input  is  used  for  recursive  identification  of  the  Preisach  density  function. 
A  practically  important  issue  is  the  choice  of  the  discretization  level  L  for  the  input.  Fig.  8  shows  the  identified 
density  distribution  for  different  discretization  levels  after  eight  periods.  The  output-based  gradient  algorithm  is 
used  with  7  =  1. 


a  (A/m) 


(a) 


(b) 


xIO4  1-5 


2 

(3  (A/m)  2  5 

3 


(C)  (d) 

Fig.  8:  The  Preisach  density  function  identified  for  different  levels  of  discretization  L.  (a)  L  =  5;  (b )  L  =  10;  (c) 
L  =  15;  (d)  L  =  25. 
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Although  it  is  expected  that  the  higher  discretization  level  L,  the  higher  model  accuracy,  there  are  two  factors 
supporting  a  moderate  value  of  L  in  practice:  the  computational  complexity  and  the  sensor  accuracy  level.  Since 
the  number  of  cells  on  a  discretization  grid  scales  as  L2,  so  is  the  computational  complexity  of  the  recursive 
identification  algorithm.  It  should  also  be  noted  that,  from  Table  1,  the  convergence  rate  y/1  —  c\  decreases  as 
L  increases.  Fig.  9  shows  the  CPU  time  used  in  recursive  identification  for  different  discretization  levels.  (To 
obtain  the  CPU  time,  the  recursive  algorithm  is  carried  out  again  using  the  collected  data  (the  current  /  and  the 
displacement  D )  of  8  periods  on  a  Dell  laptop  Inspiron  4150.)  Also  shown  in  Fig.  9  is  the  CPU  time  it  takes  to 
compute  the  Preisach  density  function  using  an  off-line,  constrained  least  squares  algorithm  [7],  where  the  data  of 
one  period  were  used.  From  Fig.  9,  the  square  law  for  the  recursive  algorithm  is  evident.  The  off-line  algorithm 
becomes  prohibitively  time-consuming  as  L  gets  large,  due  to  the  increasing  complexity  of  solving  a  constrained 
optimization  problem  of  many  variables. 

In  the  presence  of  the  sensor  noise  and  unmodeled  dynamics,  higher  discretization  level  may  not  necessarily 
lead  to  improved  performance.  Fig.  10  compares  the  measured  hysteresis  loops  against  the  predicted  loops  based 
on  the  identified  parameters  for  different  L.  Although  the  scheme  with  L  =  10  achieves  much  better  match  than 
the  scheme  with  L  =  5,  there  is  little  improvement  when  L  is  increased  to  15.  Hence  for  the  Terfenol-D  actuator 
(and  the  sensor  used),  it  is  determined  that  L  =  10  is  an  appropriate  discretization  level  for  the  Preisach  operator. 


Fig.  9:  CPU  time  used  in  the  recursive  gradient  algorithm  and  the  off-line  least  squares  algorithm. 


4  Adaptive  Inverse  Control 

Fig.  11  shows  a  schematic  of  adaptive  inverse  control.  The  input  to  the  Preisach  operator  is  obtained  through 
inversion  of  the  Preisach  operator  T  with  the  current  estimate  of  the  Preisach  density.  The  error  between  the 
reference  trajectory  and  the  achieved  trajectory  is  then  used  to  update  the  parameter  estimate  and  thus  the 
inverse  model.  For  a  Preisach  operator  with  piecewise  uniform  density,  its  exact  right  inverse  can  be  efficiently 
constructed  [13].  Note  that  for  the  Preisach  operator  T  (with  a  piecewise  uniform  density  function)  to  have  a 
unique  (right)  inverse,  it  is  required  that  all  diagonal  cells  have  strictly  positive  density  values  [17].  However, 
when  a  particular  inversion  algorithm  (e.g.,  the  one  in  [13])  is  used,  the  inverse  trajectory  can  be  made  unique 
even  if  some  diagonal  density  values  are  zero. 

To  the  authors’  best  knowledge,  the  parameter  convergence  issue  in  adaptive  inverse  control  of  hysteresis  has 
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1.5 


(a) 


(b) 


(c) 


Fig.  10:  Comparison  of  measured  hysteresis  loops  with  predicted  loops  based  on  the  identified  Preisach  density 
function,  (a)  L  =  5;  (b)  L  =  10;  (c)  L  =  15. 


not  been  studied  in  literature.  In  this  section  we  aim  to  answer  the  following  question:  how  is  the  parameter 
convergence  related  to  the  reference  trajectory? 


Fig.  11:  A  schematic  of  adaptive  inverse  control. 


Theorem  4.1  Assume  that  the  true  Preisach  operator  has  a  piecewise  uniform  density  function,  as  illustrated  in 
Fig.  3  (a).  Denote  by  v*  the  K -dimensional  vector  of  true  densities,  and  by  ysat  the  saturation  output  corresponding 
to  v* .  Let  the  output-based  gradient  algorithm  be  used  for  the  parameter  update.  Then 

(1)  For  any  reference  trajectory  yref\m]  with  yref\p\  €=  [— ysat , ysat] ,  the  parameter  estimate  v[n \  — >  f>oo  for  some 

f>oo;  and  the  tracking  error  e[n\  =  y[n\  —  yref[n \  — >  0  as  n  —>  oo; 

(2)  Assume  that  the  density  of  the  cell  (1,1)  is  positive.  Let  yref\’\  be  periodic  of  period  N  that  visits 

—Vsat,  and  without  loss  of  generality  yref[  1]  =  — ysat •  Define  ur[n]  =  T~1[yref[’],  ^o]  [n\  with  if  o  the  memory  curve 
corresponding  to  the  negative  saturation,  where  T_1  is  as  constructed  in  [13].  Then  ur[]  is  also  periodic  with  period 
N ,  and  ur[  1]  =  ^min.  Let  the  vector  of  signed  areas  of  cells  corresponding  to  ur[f]  be  Wr[-\  (which  is  also  periodic), 
and  the  null  space  of  [Wr  [1]  •  •  •  Wr  [A^]]T  be  Afr.  Then  the  parameter  estimate  v[n\  — >  J\fr  +  n*  =  {x  +  z/*  :  x  G  Mr}. 

In  particular,  if  {Wr[n\}][=1  spans  u[n\  — >  v* .  Analogous  results  hold  if  ul,l  >  0  and  yref[m\  visits  ysat- 


Proof.  (1)  Define  v[n\  =  v[n\  —  v*,  and  5[n\ 


=  v[n]T v[n\.  From  the  output-based  gradient  algorithm  (letting  7  =  1 
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without  loss  of  generality), 


Since  S[n\  >  0,  S[n\  — »  ^oo  as  n 
leads  to 


6[n  +  1] 


<%] 

5[n\ 


z>  [n] T 11^  [n]  II7  [n] T  £  [n] 


(y\n\  ~  y[n])2 
W[n]TW[n } 


<  S[n\. 


oo.  This  immediately  results  in  u[n\  — >  for  some  v^. 


Y^(y[n]-y[n])2 

^  W[n]TW[n } 


=  <S[0]  -  Soo  <  oo, 


which  implies  0.  Since  IF[n]TIF[n]  <  C  for  some  constant  C, 


(11) 


Summing  (11)  over  n 


| y[n\  —  y[n]  \  — >  0,  as  n  — >  oo. 


(12) 


The  Preisach  operator  f  based  on  parameter  estimate  at  n  can  be  exactly  inverted  (hence  y[n\  =  yref  [n])  except 
for  the  following  two  cases:  (a)  ysat[n }  <  yref[n\  <  ysat ,  and  (b)  -y8at  <  yref[n]  <  - ysat[n ],  where  ysat[n } 
denotes  the  saturation  output  corresponding  to  0[n\.  For  the  case  (a),  the  inversion  cannot  be  exact  and  the  input 
under  the  inversion  algorithm  is  u[n\  =  umax  (so  that  the  corresponding  y[n\  is  closest  to  yref[n]),  which  implies 
y[n\  =  ysat[ri\,  y[n\  =  ysat ,  and  hence  | y[n\  —  yref  [n] \  <  \y[n]  —  y\n\\.  The  same  conclusion  holds  for  the  case  (b). 
It  then  follows  from  (12)  that  the  tracking  error  e[n\  =  y[n\  —  yref\p\  approaches  0  as  n  — »  oo. 


(2)  Since  yref\kN  +\]  =  —  ysat  for  each  k  and  >  0,  ur[kN  + 1]  =  and  the  state  of  the  Preisach  operator 
is  reset  at  kN  +  1.  The  periodicity  of  ur[]  then  follows  from  that  of  yref[m]  and  the  inverse  algorithm.  From  the 
first  part  of  the  theorem,  \y[n\  —  yref  [n]\  — >  0  and  hence  y[kN  +  1]  —>  — ysat  as  k  — >  oo.  Again  from  z/yi  >  0, 
the  input  u[kN  +  1]  approaches  umin  as  k  —>  oo,  and  u[kN  +  m\  —>  ur[m\  for  1  <  m  <  N.  As  a  consequence, 
W[kN  +  m\  Wr[m\.  Since  as  k  oo 

W[kN  +  m]T v[kN  +  m\  — >  0,  for  1  <  m  <  TV, 

we  conclude  Wr  [m]Tz>oo  =  VFr[m]Tz/*,  1  <  m  <  AT,  i.e.,  z>oo  G  A/"r  +  z/* .  Analogous  arguments  can  be  used  for  the 
case  where  T>l,l  >  0  and  yref[]  visits  ysat •  □ 


Simulation  has  been  conducted  to  illustrate  Theorem  4.1.  Fig.  12  shows  the  simulation  result  of  tracking  a 
sinusoidal  signal  with  amplitude  ysat  using  the  output-based  adaptive  inverse  scheme.  One  can  see  that  the  tracking 
error  goes  to  zero.  Fig.  13  compares  the  parameter  estimate  after  160  periods  with  the  true  parameter  values. 
Note  that  the  individual  density  values  do  not  converge  (Fig.  13(a)  ).  For  this  particular  reference  trajectory,  the 
asymptotic  input  will  be  periodic  varying  between  umin  and  rzma;E  without  other  reversals.  Hence  for  each  i  the 
signed  areas  of  cells  (i,  1),  •  •  *  ,  (i,  i  —  1)  will  be  the  same;  and  for  each  j,  the  signed  areas  of  cells  (j  +  1,  j),  •  •  •  ,  (L,  j) 
will  be  the  same.  What  separates  a  diagonal  cell  from  other  cells  of  the  row  (or  the  column)  is  its  triangular  shape. 
As  a  result,  one  can  predict  that  the  densities  of  diagonal  cells  will  be  correctly  identified  and  the  sum  of  densities 
of  cells  in  each  row  (or  column)  (excluding  the  diagonal  element)  will  be  correctly  identified.  This  is  verified  in 
Fig.  13(b). 


Experimental  results  for  tracking  a  sinusoidal  signal  are  shown  in  Fig.  14  with  two  different  adaptation  constants 
7.  When  7  is  bigger,  the  trajectory  converges  to  the  steady  state  faster  but  with  larger  tracking  error  due  to  higher 
sensitivity  to  the  noise.  For  7  =  0.2,  the  tracking  error  is  under  1  yam.  Considering  the  sensor  precision,  almost 
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perfect  tracking  is  achieved  for  the  full  operational  range  of  the  actuator  (Fig.  15).  In  these  experiments  and  other 
experiments  reported  hereafter,  the  discretization  level  L  =  10. 

If  the  reference  trajectory  does  not  cover  ±ysat,  Theorem  4.1  says  that  the  tracking  error  still  goes  to  zero, 
but  one  cannot  say  more  about  the  parameter  convergence.  In  this  case  there  is  no  reset  mechanisms  during  the 
adaptation,  and  depending  on  the  initial  conditions  of  the  system  and  the  adaptation  process,  the  final  steady- state 
input  trajectories  can  be  different  (while  the  output  trajectories  are  all  consistent  with  the  reference  trajectory). 
Essentially  in  this  case  there  may  exist  multiple  minor  loops  that  satisfy  the  output  requirement.  This  is  also 
demonstrated  by  experiments  (Fig.  16),  where  two  different  inputs  are  found  to  achieve  tracking  of  a  sinusoidal 
trajectory  with  amplitude  15/im  (and  a  DC  offset  of  30/im). 

We  note  that  in  Theorem  4.1  the  output-based  algorithm  is  required.  The  results  cannot  be  extended  to  the 
difference-based  algorithm.  Indeed  it  has  been  verified  that  the  difference-based  adaptive  inverse  scheme  cannot 
achieve  asymptotic  tracking  (Fig.  17). 


5  Conclusions 


This  paper  has  been  focused  on  recursive  identification  and  adaptive  inverse  control  of  hysteresis  in  smart  materials. 
A  Preisach  operator  with  piecewise  uniform  density  function  was  used  to  approximate  smart  material  hysteresis.  On 
the  theoretical  side,  a  necessary  condition  and  a  sufficient  condition  for  the  parameter  convergence  were  presented 
in  terms  of  the  input  to  the  Preisach  operator.  In  contrast  to  the  results  for  linear  systems,  the  conditions  here 
center  around  the  local  maxima/minima  of  the  input.  Asymptotic  tracking  under  the  output-based  algorithm  was 
established,  and  the  issue  of  parameter  convergence  was  discussed  in  terms  of  the  reference  trajectory. 

Practical  issues  of  using  the  adaptive  schemes  were  studied  in  depth,  through  both  simulation  and  experiments. 
Two  types  of  adaptive  gradient  identification  algorithms  were  compared.  It  was  found  that,  for  purely  recursive 
identification  (i.e.,  given  the  same  input),  the  difference-based  method  has  a  higher  convergence  rate  in  the  absence 
of  noise,  but  it  is  more  sensitive  to  the  measurement  noise.  On  the  other  hand,  for  tracking  with  adaptive  inverse 
control,  the  difference-based  scheme  is  not  usable.  The  choice  of  the  level  of  discretization  was  also  discussed.  It 
was  shown  that  a  moderate  L  can  be  used,  taking  into  account  the  factors  of  computational  cost,  sensing  precision, 
and  model  accuracy.  In  particular,  experimental  results  have  shown  that  the  adaptive  scheme  is  able  to  virtually 
cancel  out  the  hysteresis  effect  throughout  the  actuator  working  range  with  L  =  10  for  our  experimental  setup. 
Most  results  of  this  paper  are  applicable  if  the  adaptive  least  squares  algorithms  [20,  21]  (instead  of  the  gradient 
algorithms)  are  used  for  parameter  update. 

For  future  work,  it  will  be  of  interest  to  extend  the  results  reported  here  to  the  cases  where  the  hysteresis 
output  is  not  directly  measurable.  Such  cases  happen  if,  e.g.,  the  high-frequency  dynamics  of  the  smart  material 
actuator  is  not  negligible  [13],  or  the  actuator  is  used  to  control  some  other  plant. 
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Fig.  12:  Simulation  result  of  adaptively  tracking  a  sinusoidal  signal  with  amplitude  ysat  (L  =  10). 


18 


(a)  (b) 


Fig.  13:  Comparison  of  identified  parameter  values  with  true  values  L  =  10:  (a)  Comparison  for  individual  cell 
density  values;  (b)  Comparison  for  aggregate  cell  density  values. 


(a) 


(b) 


Fig.  14:  Experimental  results  of  tracking  a  sinusoidal  reference  trajectory:  (a)  7  =  0.2;  (b)  7  =  0.5. 
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Fig.  15:  Achieved  displacement  vs.  desired  displacement  (over  the  full  operational  range  of  the  actuator). 


Fig.  16:  Two  minor  loops  corresponding  to  tracking  a  sinusoidal  signal  with  amplitude  less  than  ysat. 


Fig.  17:  Experimental  result  showing  that  the  difference-based  adaptive  scheme  cannot  achieve  asymptotic  track¬ 
ing. 
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